Exact results in the large system size limit for the dynamics of the Chemical Master Equation, a one 

dimensional chain of equations 
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We apply the Hamilton- Jacobi equation (HJE) formalism to solve the dynamics of the Chemical Master Equa- 
tion (CME). We found exact analytical expressions (in large system-size limit) for the probability distribution, 
including explicit expression for the dynamics of variance of distribution. We also give the solution for some 
simple cases of the model with time-dependent rates. We derived the results of Van Kampen method from HJE 
approach using a special ansatz. Using the Van Kampen method, we give a system of ODE to define the vari- 
ance in 2-d case. We performed numerics for the CME with stationary noise. We give analytical criteria for the 
disappearance of bi-stability in case of stationary noise in 1-d CME. 

PACS numbers: 02.50.Ey, 87.10.-e 



I. INTRODUCTION 



The statistical physics of a living cell requires a theory 
for chemical reactions with few molecules Q]],£|]. One of 
mathematical tools here is the Chemical Master Equation 
(CME) HH, describing the dynamics of probability P(X, t) 
of having different (integer) number X of molecules. The 
same CME equation could be applied to other areas of sci- 
ence as well, even in the financial market theory J15|]. Re- 
cently in |@] has been considered a kinetic model [7] for the 
phosphorylation-dephosphorylation cycle in the cell, and the 
corresponding CME was investigated. The authors considered 
the model with the bi-stability phenomenon, and have been 
derived some results for the existence of bi-stability. Here we 
solve exactly the dynamics of CME, a very important topics 
according to |g]. The accurate (exact) solution of the CME 
dynamics is important for financial market modelling @]. In 
this article we will derive exact dynamics for the CME. 

The master equation is formulated as a system of linear dif- 
ferential equations for P(X, t)of having X molecules, < 
X < N, where N is a large integer: 

d -^ = R + (^l)P(X-l,t ) + 
Ndt N ' ' 

R-(^-)P(x + i,t) - (#+(§) + R-(^))P(x,t)(i) 

Here R + is the growth rate and i?_ is the degradation rate. 
Actually we should modify the equation at the border: 

= R_(±)P(1, t) i? + (0)F(0, t) (2) 



and 

dP{N, t) 



XiU =R + (l-±)P(l-±t)-R-(l)P(N,t) (3) 



The large parameter N describes the system volume. 

A close master equations have been considered in evolution 
theory in H-H. 

Let us introduce the coordinate x and function p(x, t): 

x = X/N, p(x, t; N) = NP(X, t) (4) 

Assuming that the probability distribution is a smooth func- 
tion of x, 



P(X + l,t)-P(X±l,t)<^ 1 



(5) 



one gets the Fokker-Plank equation for the model by Eq.(l). 

The investigation of CME via Fokker-Plank equations meet 
some problems [8]. An alternative approach is to assume that 
p(x,t) is not a smooth function of x, but the function u(x, t) 
is, where 



p(x,t;N) = exp[Nu(x,t)] 



(6) 



'Electronic address: saakian@yerphi.am 



Thus the p(x, t; N) might be un-smooth in the limit of N —> 
oo, but still have smooth u(x,t). Such an ansatz has been as- 
sumed for the statics in fl, and for the dynamics in lfl2[[l3ll . 
while considering the evolution models. This ansatz gives the 
solution of the dynamics and steady state with the accuracy 
0(1/N), while the approximation of the master equation via 
Fokker-Plank equation, assuming the smoothness of p(x, t), 
is certainly wrong. This topic we already discussed in IU4I1 . 
then in ITai . 

In section II we calculate the dynamics of variance for pop- 
ulation distribution, using HJE. In section III we solve the 
CME dynamics for the simple case of time dependent rates. 
In appendix we re-derive the variance of distribution using 
Van Kampen method, also used that method to calculate the 
variance for general 1-d multi-step CME, as well as give ODE 
to derive the variance in case of 2-d CME. 
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II. THE MASTER EQUATIONS WITH CONSTANT RATES 

A. Hamilton-Jacobi equation for the Chemical Master 
Equation 

Using ansatz by Eq. (6), the model equations (1) can be 
written as Hamilton-Jacobi equations for u = \np(x, t)/N 



Fin 

^+H(u',x)=0 > 



(7) 



where vl = du/dx, 

H(u', x) = R+ 0) + R- 0) - R+ {x) e - u ' - i?_ (x)e u ' (8) 

Eqs. (7), (8) has been derived in lfl6ll . lfl8ll . where has been 
investigated mainly the corresponding Hamilton equations to 
calculate the time of extinction from meta-stabile state. We 
are interested in the investigation of the whole distributions, 
using the traditional mathematical method of characteristics. 
To solve the HJE, we consider the Hamilton equation for x and 
corresponding momentum, getting the system for the charac- 
teristics lfl9ll20ll 

x = H v [x,v) = R+(x)e- v - R_[x)e v , 
v = -H x (x,v), 
it = v H v (x,v) — H(x,v) — vx + q, (9) 

subject to initial conditions: ir(0) = xo, v(0) — vq(xo), 
u(0) = uq(xq). Here v := du/dx, vq(x) := u' (x), q := 
du/dt. 

The respective solution to Eq. (0 in (x, t) - space is called 
the characteristic of Eq. (Q. 

Our Hamiltonian is time independent. Then Eq. (0 and Eq. 
(O result in 







(10) 



Along the characteristic x = x(t) the variable q is constant, 
so q is selected to parameterize these curves. 
Consider the equation 

q = -(R+(x) + R-(x)) + R+{x)e- v + R_(x)e v (11) 
It has a solution for 

q > 0, (12) 

if at some point 

R+(x) = R.(x), (13) 

and we take q > 0. 

Using Eq. (1 1), we transform the first equation in (0 into 



x = ±y/(q + R+ + R^f - m+R. 
Consider the following initial distribution 

u${x) — —a(x — xq) 2 . 



(14) 



(15) 



with large a. 

The maximum of the distribution corresponds to the point 
u' = 0, therefore q = 0. 

Thus for the maximum of distribution we should consider a 
characteristic with q = 0. 

Integrating the Eq.(14), we derive: 



T = 



dy 



R+(y)-R.(y) 



(16) 



Such equation was derived in |g]. 

We can define the dynamics of the full distribution. 
Let us define the function: 



T(x,q) 



dy 



,o V(<i + R+(v) + R-(vW - 4R+(y)R-(v) 

(17) 

To calculate u(x, t) we first calculate q for the given x from 
the equation 

T{x,q)=t (18) 
Eq.(18) defines an implicit function 

q=Q(x,t) (19) 

B. The elasticity 

It is important to calculate u"(x, t) at the point of the max- 
imum of distribution, the "elasticity" 

We calculate q' x = d 2 u/dxdt from Eq.(19): 



dq dQ{x,t) 
dx dx 



(20) 



We can calculate the last derivative at fixed t from the expres- 
sion: 



T(x, q) = const 
dT{x,q) dT{x,q) , 

dx dq q ' x 



(21) 



It is equivalent to calculate q' x from Eq. (17) at the fixed t. 
Thus we get the following equation for q' at the point of 
maximum (u'(x, t) = 0): 

<f r ^#4^8 = 0(22) 



R + (x)-R-(x) - 1 J X0 [(R + (y) - R-(yW 

From Eq.(l 1) we can obtain: 

q ' x = -(R + (x) - R-(x))v' (23) 

Thus eventually we get: 

-1 , H M D * „2 f x dy(R + (y) + R-(y)) .... 
— = {R+{x) - R-[x)) / (24) 

Eq. (24) is the main result of our work. 

In Fig. 1 is given the comparison of analytical result with 
the numerics. 
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C. Probability distribution 



Eq.(17) defines the q(x,t) for given x,t, we can define 
v(x, t) also using Eq.(7). 

To calculate u(x, t) at the given (x, t), let us consider the 
trajectory of points x(t),t) connecting that point with the 
starting point (xq, 0). 

We have chosen the trajectory to have g(a;(r),r) = q. We 
take x(t) just as a solution of the equation 



T = T(x(r),q) 



(25) 



At any point of our trajectory x(t),t), we can calculate 
v(x(t), t), while q(x(T), t) = q is constant. Eq.(ll) gives 

R_e 2v ~(q + R + + R_)e v + R + = 0, 
v(y)=- In 2i?_ + ln[((g + R+ (y) + i?_ (y) ± 

y/(q + R+ (y) + R- (y)) 2 - 4R+(y)R- (y))} (26) 

where we denoted x(t) = y. 

We derive the solution of the original Eq. (0, integrating 
the equation u — vx + q along our trajectory (the characteris- 
tics connecting the points (x. t) and (xq, 0)): 



u(x, t) = u{x, 0) + / dyv(y) + qt = 

J xo 

f dy[~ In 2i?_ (y) + ln[((q + R+ (y) + R_ (y) ± 

V(<z + R+(y) + R-{y)) 2 - 4R+(y)R-(y))}} + qt (27) 

Having the expression u(x, t) we can calculate p(x, t). 



D. The restricted meaning of probability distributions in 
master equation 

In case of evolution models [9-13], we have master equa- 
tion similar to Eqs.(l)-(3), only the negative term ~ P(X, t) 
has other coefficient than — + i?_ ), and therefore there is 
no a balance condition 

Contrary to the case of master equations in evolution mod- 
els j^- lflill where all the initial distributions have a mean- 
ing, now there are some restrictions. We should clarify well 
the meaning of the probabilities P(X,t). At every mo- 
ments of time the system has only ONE value of X, and 
P(X, t) just gives such probabilities. We should solve the 
system (1),(2),(3) for the given initial value P(Xo, t) = 1 and 
P(X, t) — for other X, other initial distributions have not 
meaning. 

Another difference is connected with the stable point solu- 
tions. Eq. (13) gives the steady state solution. If that equation 
has two stable solutions x i,X2 and the probability of these 
two positions is the same, i.e. 



* 2 , . R+(x) n 
X1 R-{x) 



(28) 



(t) 



for the model 



FIG. 1: The graphics for the elasticity V(t) 
withiV = 100,R+(x) = exp(-x),R+(x) = exp(x),x{0) = 0.5. 
The smooth line is the analytical result by Eq.(24) and the dashed 
line is the numerical result. The difference is less than 0.5%. 



at fci = 43.1274. One can use the HJE method to calculate 
the mean period of time that solution, initially located at one 
stationary point, will move to the other stationary point lfl7ll . 
Then it again should return back, as every moment the system 
could exists only with one value of X. 

In case of evolution model, the system goes to the equilib- 
rium state instead of oscillation between two stable solutions. 



E. The dynamics for the stationary but random rates. 

Consider now the case when the rates are smooth functions 
of x plus some random noises. The noise in the rates is well 
confirmed experimentally [2j]]. We took the case of rates from 

a 

R + (k u x) = (1 - x)(0.5 + k lX 2 ), 

R-(k 2 ,x) = x(k 2 +O.OI2 2 ) (29) 
where now k±9 and k 2 are random variables, 

ki = Kiexp(a^i(x) - a 2 /2)x 2 ), 

k 2 = K 2 exp« 2 (x) -a 2 /2) (30) 

^i(x),£, 2 {x) have normal distributions. We consider the 
model with N = 100, and performed numerics for different 
values of parameters K\, K 2 ,a. For K\ = 50, K 2 = 10, at 
a k, 0.9 there is a phase transition: instead of two steady state 
solution we get one steady state x ~ 0.088. 

We can analytically estimate the transition point (from one 
stabile point to bi-stability), considering the behavior of 



U{x) =< In—— > | fel)fe2 



(31) 



we again should accurately interpret the HJE results @]. For 
the rates given by Eq.(27) at a — 0, k 2 — 10, the transition is 



R-{k 2 , x 

We found that the function U(x) changes its behavior with the 
level of the noise. At a = it has only three roots U (x) = 0, 
while for a > 0.75 there is one root. 



m. MASTER EQUATION WHEN THE RATES VARY 
WITH TIME 



A. The simplest solvable case 

Consider the case when birth and death rate coefficients in 
CME change with the time as g(t) and f(t), while they are 
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U(x) 



0.2 0.4 0.6 XQ.8 1.0 



It is interesting to find the variance of distribution. Differenti- 
ating Eq. (36), we get, putting vq = at the maximum point: 

- v'(x, t) = — f dr (39) 

IoW)+f(r)] 

Thus the variance of distribution decreases with the time. 



FIG. 2: The graphics for the U(x) = ]n(R+(x)/R-(x)) at K 1 = 
44, Ka = 10, a — 0. There are three solutions for the equation 
U(x) = 0. 



U(x) 



fl:2 0.4 0.6 X0.8 1.0 



FIG. 3: The graphics for the U(x) =< la(R + (x)/R-(x)) > at 
K\ — 44, K2 = 10, a — 0.8. There is only one solution for the 
equation U(x) = 0. 



the same for all the x. We have the following Hamiltonian: 

- H = g(t)e~ v + f(t)e v — g(t) — f(t), v = — (32) 

We get for the characteristics the following system of equa- 
tions: 



dp dH 



= 



(33) 



dt dx 
Thus the p is constants: 

v = v , (34) 

For the coordinate we have a simple equation: 

dx _ dH 
dt dv 
Thus we have a simple solution 



g(t)e- v ° - f(ty 



(35) 



x-x = e~ v ° I g(r) - e Va / f(r) (36) 



n 



(i 



We find vo as a function of x, t from the solution of the last 
equation. 

u(x,t) = u(xq,0) + v (x — x ) + 



[g(r)e- V0 + f(r)e va - 9 (t) - f( T )]dr (37) 



For the maximum of distribution we have vo = 0, therefore 
we get an equation: 



x-x = / [g(r) - f(T)]dT 
Jo 



(38) 



B. Rates as linear functions of the number of molecules 

Consider now the Hamiltonian 

-H = (ax+g(t))e- v + (bx + f(t)y-(a+b)x~f(t)-g(t) 

(40) 

We have 

'''' ae~ v + be" - [a + b) 



dt 

We can solve this case: 



(41) 



dv 



ae~ v + be" - (a + b) 

(42) 

We can define the function v = V(vq, t) from the latter equa- 
tion 

~ae at - ae bt - be at+vQ + ae bt+v0 



V(v ,t) = Log 

We have also 
dV 



dv t 



■(*) 



ae at _ b( ,bt _ 5 e at+v() + fr e fct+v0 
_ be at+v + ae bt+v 



(43) 







ae at - ae bt - be at+v ° + ae bt+v ° 



-be at+v ° + be bt+v ° 



(44) 



(45) 



ae at - be bt - be at+v ° + be bt+v ° 
Now we solve the equation for x: 

^ = (b-a)x + f(t)e v -e- v g(t) 
Its solution gives 

x-x = e {b - a)t f [f{T)e v ^ ^-g(T)e- v ^ V0 ^]e^ b ~ a ^ T dT 
Jo 

(46) 

Eqs.(42),(46) together define the trajectory of characteristics. 

To get the solution for the maximum, we put v = on the 
left hand side of Eq. (42). In this way we get a simple explicit 
equation for vq as a function of time t. 



ae 



at - ae bt - be at+vQ 



ae 



bt+vO 



1 



(47) 

ae at _ be bt „ b( ,at+v0 + 5gbt+v() 

Putting that solution into Eq.(46), we find the trajectory of the 
maximum of distribution. 

To calculate v' we should differentiate the expression in 
Eq.(46) via x. Using the relation V'(vo, t) = 1, we derive 

1 



e (6-«)t x 



[f(r)e 



V(v Q , 



,dV 



■ g ( T )e- vivo ' r) ] — (r)e^ (b - a)r dr (48) 
dv 



where vq is defined by Eq.(46) as a function of t. 



5 



IV. CME IN MULTI-DIMENSIONAL SPACE 

Consider now the CME in multi-dimensional space, when 
we have d-dimensional X. 

dP(X,t) 



K 



]T Rn£^)P{X n,t) - £ifc(£)P(X,t)(49) 



-K 



Ndt 

I' 



We get HJE for the function u(x) with following Hamiltonian: 

dH 



dt 



+ H(x,p)=0, 



K 



-H(x,p)= J2 ^ s (f)[exp[-np]-l] (50) 
m=-K,l<l<d 

Let us investigate the motion of the maximum of distribution. 
We denote y(t) = x the maximum of distribution at moment 
of time t, and assume the following ansatz for the u(x, t) near 
the maximum of distribution: 

u{S,t) = ~ lt;x-y(t)\V\x-y(t) gt; (51) 

Putting this ansatz into 

< t + H' xi (x,p) + Y,H' f Jx,p)^ = (52) 

m 

Consider the point x = y(t),p = 0, and use the identity: 

dp? 



We get: 



Y,V lm ^-( Rm...n d J2n m )V lm =0 (54) 



dxi 



K 



(53) 



dt 



n h = -K 



We get the following ODE for the dynamics of the maximum 
of distribution: 

—fa- = Z^Qmy m {t)), 



l<h<d -K<n h <K 



V. DERIVATION OF ELASTICITY IN 1-D CASE 

Let us get system of ODE for the V. We consider the multi- 
step version of CME in 1-d. 

dP(X, t) _ 
Ndt 

Rn£^)P{X - n ,t)-J2R n d)P(X,t)(56) 



n=-K 



Now we have the following Hamilton-Jacobi equation: 

K 

-H(u',x)= Rn{x)[e- nu ' -1] (57) 

n=-K 

Now we have to consider the higher terms in expansion of 
u(x, t) near the maximum. 

u = -V(x - y(t)f/2 - T(x - y(t)) 3 /6 (58) 

Eq.(55) gives 

dy(t) 



dt 



y(t))b 



b= Y R n{x)l 



-K<n<K 



We also denote: 



a{x) = Y R n(x)r 



(59) 



(60) 



—K<n<K 

therefore 

y = b 

With the ansatz by Eq.(56) we have: 

< xt = ~V + Tb, 
u^ t = 2Vb-Tb 2 + Vy = 
2Vb - Tb 2 + Vbb' 



(61) 



From the other hand we have differentiating the right hand 
sideofEq.(50): 

- uZ t = H>; p V 2 - 2H'J p V - TH' p - aV 2 - 2b'V - bT, 
-uZ t = -K p Vy + H' p ' p bV 2 - H' p V + H p Tb 2 = 

-bb'V + abV 2 -bV + Tb&2) 

We derived Eqs.(58),(59) putting x = y{t). 
Then 

- V + 2Tb = aV 2 - 2b'V, 
3Vb - 2Tb 2 = aV 2 b - 2b'V 



(55) Removing T we get 



or 



2b— V = -Abb'V + 2baV 2 
dt 



d 1 ,/ 1 „ 

dtv= b v- 2a 



(63) 



(64) 



Which gives Eq.(24) for 1-step CMW case. 

For the multi-step CME we have the following expression 
for the elasticity: 



dyajy) 
[b(y)} 3 



(65) 
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VI. CONCLUSION 

In conclusion, we calculated exact probability dynamics 
Eq.(27) for the Chemical Master Distribution using Hamilton 
Jacobi equation method. Using the methods of characteristcs 
for solution of HJE, we gave explicit expression for the vari- 
ance of distribution Eq.(24). The latter is important both for 
chemical yO and financial applications. We derived the vari- 
ance also directly from HJE, using a special ansatz. The latter 
is equivalent to Van Kampen method 04J]. Using Van Kam- 
pen method, in Appendix we give an exact expression for the 
variance for general 1-d model, as well as a system of ODE 
to define the variance in 2-d case. Both HJE and Van Kampen 
methods gave identical results for the variance and the maxi- 
mum point dynamics of CME: HJE gives directly differential 
equation for the vasrianc, while the Van Kampen method orig- 
inally gives a differential equation for the probability distri- 
butions, Eq.(A.5), which later proved the differential equation 
for the variance. HJE gives also exact steady state distribution, 
and is more adequate for investigation of meta-stabile points 
fl^l Jl8ll . Using HJE, we derived exact dynamics Eqs.(44)- 
(46) for the case of 1 -d CME when the rates are linear func- 
tions of the coordinate (number of molecules) plus some time 
dependent functions. 

We performed some numerics in case of static noise, also 
give an analytical criteria for the level of the noise when the bi- 
stability disappears. Our choice of potential for the averaging 
Eq.(31) is rather arbitrary (contrary to all other results in this 
article, which are derived rigorously and are exact). A further 
investigation of the problem is necessary. Perhaps it is pos- 
sible to investigate the more realistic case of non-stationary 
noise 112111 . 

DBS thanks DARPA Prophecy Program and Academia 
Sinica for financial support. 



Appendix A: The calculation of variance dynamics using the 
Van Kampen method 

1. 1-d multi-step models. 

Consider the CME defined trough the system of equations: 

dP(X,t) _ 



K 



£ R l (^±l)P(X + l,t) -£i*(^)P(X,i) (A.1) 



l=-K 



Ndt 
X 

* TV- 



Following to Van Kampen, we assume that near the maximum 



(A.2) 



X = N<f>(t) + VN£{t), 

p(x,t) = n^,t) 



For the maximum point </>(i) = ^ iP(i, t) we will derive an 
equation (4) later. 



We consider terms with different scaling by N in the equa- 
tion: 



K 



n y, 



N N 



)P(X-l,t) 



-Nj2MHt) + ^)P(X,t) (A3) 



— K 



Collecting together the ~ V A*" terms, we get 

d(f>(t) h 



dt 



= ]T mm)) 



(A.4) 



l = -K 

The N° terms give an equation 

dU(t) v a d 2 _ 

b = Y lR h a = ^Z l2R i 

i i 

We derive the following equations for < £j >: 

d<£> 



dt 



b'(t) < & > 



(A.5) 



(A.6) 



The initial condition < £(0) >= gives < >= 0. 
We get for the variance 

j t < m 2 >= <m) + xm) < m >? (a.?) 

We can consider ODE via ip instead of t. Then Eq. (A.4) gives 
dip/dt = b, then: 

b-^- < e >= a(V) + 2&'(V) < £ >) 2 (A.8) 
dtp 

Eventually we get Eq.(24), with R + — i?_ — > b = 
Y^i Rd- R+ + R- a = Y^,i I Ri- 



2. 2-d case 

We will apply the Van Kampen method, giving ODE to cal- 
culate the variance of distribution in 2-d case. 

Consider the following system of equations for 

P(X,Y,t),0 < X < N,0 < Y < N: 

dP(X,Y,t) X Y X X 
dt - _ A' N> + Rl - { N ] > N ] 

RM^-^)P(x~W) + 

Rr-(^,^j^)P(X + l,Y,t) 
R2 + (§,^-)P(X,Y-l,t) + 

R 2 -(§,^P : )P(X,Y + I,t) (A.9) 
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Following to [4], we introduce the fluctuating variables £1, £2 
and replace P(X, Y, t) by nfo, £ 2 > *),see ff: 



X = JV0 1 (t)+VJVei(t) J 
Y = N^{t) + VN&{t), 
P{X,Y,t)=Il(£ 1 ,&,t) (A. 10) 

where ip\ (t) , (t) give the solution in case of infinite N: 

#i(t) 



dt 
dMt) 
dt 



b 1 (Mt),Mt)) 



62(^1,^2) = R2+(i>l,^2) - R2-(lpl,1p2) 



(A.ll) 



We can solve Eq.(A.l 1) and calculate ip\ (t),ip2(t). 

Following to the methods of [4J], we derive an equation: 



czn(f) 

dt 



i2 / a-)o7-(ean(Cl,6,t)) + 
Ot, a 



2 

db a d 1 ^ 9 2 



n (A. 12) 



where we denoted R' a± = dR ^^ 2{t) \ R a± 
We derive the following equations for < ^ >: 



(A. 13) 



The initial condition < ^(0) >= gives < >= 0. 
We get for the variance 



dt 



< m 2 >- 



a l (Mt)^2(t)) + 2b' l (Mt),Mt)) < >) 2 (A.14) 



We can calculate the variance numerically as function of t, 
using the solution of Eq.(A. 11). 
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